/*  sdfexpbs.cpp  by K.Tsuru  */
// function ID = 3318 ver. 2.18
// First version has been finished on 30 May 2003.
// rewrite on 4 Aug 2016 since version 2.30
/***************************************************
It evaluates exp(x) using the binary splitting method.
****************************************************/

#ifndef SN_H
#include "sn.h"
#endif

//static const SDouble ONE(1.0);
static const SLong S_BASE = Lpow10(4*DFIGURES); // consider efficiency of FFT mult.

SDouble ExpBS(const SDouble& x) {
// Check the sign of x.
    if( !x.Sign(3318) ) return ONE; // x = 0
    if( x.IsMLT(ONE) ) return (ONE + x);    // |x| << 1.0, x*x = 0;
    SDouble X(x);
    bool posX = true;
    if(X.Sign() < 0) {
        posX = false;   X.ChangeSign(3318);
    }
// Here X > 0.
    ExpArgCheck(X, 3318);   // over/underflow check ver 2.18

    SNBlock <SLFraction> slr;
    int k = SplitSDouble(X, slr, S_BASE);   // very fast
    long ip = (long)slr(0).num.Low2();

    SDouble r = Expower(ip);

    SLong num, den, a, c;
    int i = 1;

/***********************************************
The first factor is evaluated by
 exp(num/den) = {exp(num/(den*divisor)}^divisor .
************************************************/
    static const int pow2 = 16;
    static const long divisor = 1L << pow2;
    num = slr(i).num;   den = slr(i).den * divisor;
    i++;

    ExpBSR(num, den, a, c);  // exp(num/den) = a/c (a rational number)

    SDouble f = SDouble(a)/SDouble(c);  // = a/c
    f = Dpow(f, divisor);
    r *= f;

/********************************************
My note : Do not apply the fixed point mode,
 because summing up is not used,
SDouble version is faster than SLong one.
*********************************************/
    SDouble rA(ONE), rC(ONE);

    for( ; i < k; i++) {
        num = slr(i).num;   den = slr(i).den;
        ExpBSR(num, den, a, c);
        rA *= SDouble(a);   rC *= SDouble(c);
    }
    r *= rA;
    if(posX) { // r*(rA/rC)
        return r / rC;
    } else { // rC/(r*rA)
        return rC / r;
    }
}